function sineswp
%   the sineswp function computes the filtered peak amplitude based on swept sine 
%   time history data.
%
%Author:    Donald J. Hershfeld
%   ManTech International Corporation
%Created:   October 8, 2002
%Released:
%Purpose: The sineswp function requests the user to select a channel of data to be used
%         as the frequency reference. The frequency reference channel is analyzed using a Hilbert
%         transform to determine the rate of phase angle change with respect to time. 
%         The user is then prompted to select break points on the phase angle plot.
%         Logarithmically swept cosine and sine references are then fitted between the
%         break points. 
%
%         By multiplying each time history by the cosine and sine references, exactly half
%         of the energy at the reference frequency is transposed to zero frequency. The other half of
%         the energy is transposed to twice the reference frequency. A low pass filter is
%         used to create a tracking filter. Two times the output of the filter is equal to
%         the peak amplitude of the time history of the reference frequency. The cosine part
%         is the in-phase component and the sine part is the quadrature component. The cosine  
%         and sine parts are combined as a complex array in the Analysis structure.
%
%
%Function Parameters:
%   At the present time, the sineswp function does not require input parameters. However,
%   the sineswp function could be changed at a later time to require a freqrefchan.
%
%   freqrefchan:  
%       Channel number of data defined by Project.Test_Item.Run.Channel(ChanNum)to be used to
%       estimate the instantaneous rate of change of phase angle (reference frequency).
%
global Project;
% Added 29-Jul-2005 by C. McLeod
DecimateFactor = 100;
%   Request the user to select the freqeuncy reference channel.
refch=sineswpdialog;

%   Get the time reference.
TimeRef=Project.Test_Item.Run.Time_Ref;

%   Get the freqeuncy reference data.
% Modified 14-Apr-2004 by C. McLeod: Convert single to double precision
x = double(Project.Test_Item.Run.Channel(refch).Time_His);
h = hilbert(x);  %  Compute the Hilbert transform.
PhaseAngle = atan2(imag(h),real(h));  %  Calculate the instantaneous phase angle.

%  Calculate the rate of change of phase angle.
SampleRate = Project.Test_Item.Run.Sample_Rate;
freq = SampleRate * diff(unwrap(PhaseAngle)) / (2*pi);

[B,A] = butter(4, 8/Project.Test_Item.Run.Sample_Rate);
filtfreq = filter(B,A,freq);
% Added 03-May-2007 CMcLeod to eliminate values < 0 in semilogy()
filtfreqP = filtfreq;
for ii = 1:size(filtfreq,2)
    if (filtfreq(ii) < 0)
        filtfreqP(ii) = 0;
    else
        filtfreqP(ii) = filtfreq(ii);
    end
end
semilogy(TimeRef(2:length(TimeRef)),filtfreqP),grid on;  %   Display the instantaneous frequency.

set(gca,'ylim',[1 3000]);

if (Project.Test_Item.Run.Time_Ref(length(TimeRef)) < 180)
    set(gca,'XTick',0:5:Project.Test_Item.Run.Time_Ref(length(TimeRef)));
else
    set(gca,'XTick',0:10:Project.Test_Item.Run.Time_Ref(length(TimeRef)));
end

xlabel('Time - seconds');ylabel('Frequency - Hz.');

% [start stop]=start_stop_dialog;
status = 0;
while status == 0  %  Endpts not acceptable
    [xi,yi]=ginput(2);
	x = [min(xi) max(xi)];
	haxes = get(gca,'Children');
	xd = get(haxes,'xdata');
	yd = get(haxes,'ydata');
	xlow = find(x(1)>xd);
	xhigh = find(x(2)<xd);
	starttime = xd(size(xlow,2));
	startfreq = yd(size(xlow,2));
	stoptime = xd(xhigh(1));
	stopfreq = yd(xhigh(1));
	[status] = ConfirmEndpts(starttime,startfreq,stoptime,stopfreq);
end  % Endpts acceptable
first=round(SampleRate*starttime);last=round(SampleRate*stoptime);

%Modified 01-Dec-2004 by Don Hershfeld
RefFreq = filtfreq(first:DecimateFactor:last);     %  Decimate the reference frequency.

%  Calculate the referrence cosine and sine signals.
cx = 2 * cos(PhaseAngle);
sx = 2 * sin(PhaseAngle);

%  The low-pass filter.
[B,A] = butter(2, 8/Project.Test_Item.Run.Sample_Rate);

w = waitbar(0,'Processing Sine Data:  Please wait...');

% Determine Which Analysis index to which Analysis data will be saved.
try
 	n = length(Project.Test_Item.Run.Channel(1).Analysis);
 	n = n + 1;
catch
	n = 1;
end

%   Filter each channel after frequency translation.
for ich=1:Project.Test_Item.Run.Num_Chan
%     str0=num2str(ich);
%     str1=Project.Test_Item.Run.Channel(ich).Trans_Loc;
%     str2=Project.Test_Item.Run.Channel(ich).Coordinate;
%     str3=Project.Test_Item.Run.Channel(ich).Direction;
%     str4=Project.Test_Item.Run.Channel(ich).Orientation;
    % Modified 14-Apr-2004 by C McLeod: Convert single to double precision
	yr=filtfilt(B,A, Project.Test_Item.Run.Channel(ich).Time_His(first:last) .* cx(first:last));
	yi=filtfilt(B,A, Project.Test_Item.Run.Channel(ich).Time_His(first:last) .* sx(first:last));
    % Modified 03-May-2005 by Chris McLeod
	yr_mod = yr(1:DecimateFactor:length(yr));
    yi_mod = yi(1:DecimateFactor:length(yi));
	Project.Test_Item.Run.Channel(ich).Analysis(n).Type = 4; 
	Project.Test_Item.Run.Channel(ich).Analysis(n).Abs= RefFreq;
	Project.Test_Item.Run.Channel(ich).Analysis(n).Ord= yr_mod + i*yi_mod;
	Project.Test_Item.Run.Channel(ich).Analysis(n).FPoint = first;
	Project.Test_Item.Run.Channel(ich).Analysis(n).LPoint = last;
    waitbar(ich/Project.Test_Item.Run.Num_Chan);
end; % for ich=1:Project.Test_Item.Run.Num_Chan
close(w);
